function gesol = SolveFixedPoint_Model_20220926(parms,ERRTOL)
    
    MAXITERS = 10000;
    ITER_BANDS = [0;1000;5000;MAXITERS];
    UPDATEWT_BANDS = [0.9;0.95;0.98];

    NN = parms.NN;
    Nx = parms.Nx;
    Nz = parms.Nz;

    % fixed objects
    StateTransitionProbsTranspose = parms.StateTransitionProbs';
    f_fun = @(x) (1+x.^(-parms.iota)).^(-1/parms.iota);
    g_fun = @(x) (1+x.^(parms.iota)).^(-1/parms.iota);
    finv_fun = @(f) (f.^(-parms.iota) - 1).^(-1/parms.iota);
    ginv_fun = @(g) (g.^(-parms.iota) - 1).^(1/parms.iota); % Solves for theta given g
    zeta_fun = @(N,phi)(N + (phi.^(1-1/parms.chi)).*(1-N)).*((phi.*(1-N) + N).^(parms.gamma-1));
    phi_of_x = @(x)1./(1 + exp(-x));
    
    N_Nxz = repmat(parms.grid_N(:),[1 Nx Nz]);
    U_Nxz = 1 - N_Nxz;
    EXPZ_Nxz = permute(repmat(exp(parms.grid_z(:)),[1 NN Nx]),[2 3 1]);
    Y_Nxz = EXPZ_Nxz.*N_Nxz;
    sep_Nxz = permute(repmat(parms.s(:),[1 NN Nx]),[2 3 1]);
    survive_Nxz = 1 - sep_Nxz;
    x_Nxz = permute(repmat(parms.grid_x(:),[1 NN Nz]),[2 1 3]);
    phi_Nxz = permute(repmat(phi_of_x(parms.grid_x(:)),[1 NN Nz]),[2 1 3]);
    zeta_Nxz = zeta_fun(N_Nxz,phi_Nxz);
    sigma_x_Nxz = permute(repmat(parms.sigma_x,[1 1 Nx]),[1 3 2]);
    
    [N_meshgrid,x_meshgrid] = meshgrid(parms.grid_N(:),parms.grid_x(:)); % needed for interp2
    
    % dynamics for phi' = phi'(N,phi,z,z')
    x_next_Nxzz = zeros(NN,Nx,Nz,Nz);
    for iz_now = 1:Nz
        expected_znext = parms.StateTransitionProbs(iz_now,:)*parms.grid_z(:);
        if parms.OPTION_x_standardize_shock
            std_znext = sqrt(parms.StateTransitionProbs(iz_now,:)*(parms.grid_z(:).^2) - expected_znext^2);
        else
            std_znext = 1;
        end
        for iz_next = 1:Nz
            x_next_Nxzz(:,:,iz_now,iz_next) = (1-parms.rho_x)*parms.x_bar + parms.rho_x*x_Nxz(:,:,iz_now) ...
                + sigma_x_Nxz(:,:,iz_now)*(parms.grid_z(iz_next) - expected_znext)/std_znext;
        end
    end
    x_next_Nxzz = reshape(x_next_Nxzz,[NN*Nx*Nz,Nz]); % ((N,lambda,z),z') format
    phi_next_Nxzz = phi_of_x(x_next_Nxzz);
    x_next_Nxzz(x_next_Nxzz<parms.grid_x(1)) = parms.grid_x(1);
    x_next_Nxzz(x_next_Nxzz>parms.grid_x(end)) = parms.grid_x(end);
    
    % constraint on matching function (slow recovery)
    if ~isnan(parms.phi_match)
        FMAX_Nxz = 1 - parms.phi_match - sep_Nxz + sep_Nxz./U_Nxz;
        FMAX_Nxz(FMAX_Nxz>1) = 1;
        FMAX_Nxz(FMAX_Nxz<0) = 0;
        THETAMAX1_Nxz = finv_fun(FMAX_Nxz);
        GMIN1_Nxz = g_fun(THETAMAX1_Nxz);
    end
    
    % numerical constants to ensure positive consumption during
    % intermediate loops
    VMAX2_Nxz = 0.999*(Y_Nxz./parms.kappa); % C=Y-kappa*V>=0.001*Y
    THETAMAX2_Nxz = VMAX2_Nxz./U_Nxz;
    GMIN2_Nxz = g_fun(THETAMAX2_Nxz);
    
    if isnan(parms.phi_match)
        GMIN_Nxz = GMIN2_Nxz;
    else
        GMIN_Nxz = max(GMIN1_Nxz,GMIN2_Nxz);
    end

    % reserve memory for policies and value functions: ordering is (N,L,z)
    CV_S_old = zeros(NN,Nx,Nz);
    CV_S = zeros(NN,Nx,Nz);
    S_Nxz = zeros(NN,Nx,Nz);
    f_Nxz = zeros(NN,Nx,Nz);
    g_Nxz = zeros(NN,Nx,Nz);
    theta_Nxz = zeros(NN,Nx,Nz);
    V_Nxz = zeros(NN,Nx,Nz);
    C_Nxz = zeros(NN,Nx,Nz);
    
    % reserve memory for interpolation
    C_next = zeros(NN*Nx*Nz,Nz); % C'=C'((N,L,z),z')
    S_next = zeros(NN*Nx*Nz,Nz); % S'=S'((N,L,z),z')
    zeta_next = zeros(NN*Nx*Nz,Nz);
    Nnext_Nxz = zeros(NN,Nx,Nz); % N'=N'(N,L,z)
    Nnext_Nxz_notrunc = zeros(NN,Nx,Nz); % N'=N'(N,L,z)
    
    % unclassified...
    SDF_Nxzz = zeros(NN*Nx*Nz,Nz); % ((N,L,z),z') format
    StateTransitionProbs_Nxzz = reshape(permute(repmat(parms.StateTransitionProbs,[1 1 NN Nx]),[3 4 1 2]),...
        [NN*Nx*Nz, Nz]);

    % run VFI
    bContinue = true;
    N_LB = parms.grid_N(1);
    N_UB = parms.grid_N(end);

    iter = 0;
    CV_S_old(:,:,:) = reshape(reshape(max(EXPZ_Nxz - parms.b,0)...
        ./(1-parms.beta*survive_Nxz),[NN*Nx, Nz])*StateTransitionProbsTranspose,[NN Nx Nz]);

    while bContinue
        
        iter = iter + 1;

        % find policies
        g_Nxz(:,:,:) = max(parms.kappa./((1-parms.eta)*CV_S_old), GMIN_Nxz);
        g_Nxz(g_Nxz>1) = 1;
        theta_Nxz(:,:,:) = ginv_fun(g_Nxz);
        f_Nxz(:,:,:) = f_fun(theta_Nxz);
        V_Nxz(:,:,:) = theta_Nxz.*U_Nxz;
        C_Nxz(:,:,:) = Y_Nxz - parms.kappa.*V_Nxz;
        S_Nxz(:,:,:) = max(EXPZ_Nxz - parms.b + (survive_Nxz - parms.eta*f_Nxz).*CV_S_old,0);
            
        Nnext_Nxz(:,:,:) = survive_Nxz.*N_Nxz + f_Nxz.*U_Nxz; % note N'=N'(N,H,z)
        Nnext_Nxz_notrunc(:,:,:) = Nnext_Nxz;
        Nnext_Nxz(Nnext_Nxz<N_LB) = N_LB;
        Nnext_Nxz(Nnext_Nxz>N_UB) = N_UB;
        
        % interpolate: transpose data input to account for difference
        % between meshgrid and ndgrid
        for iz_next = 1:parms.Nz
            C_next(:,iz_next) = interp2(N_meshgrid,x_meshgrid,C_Nxz(:,:,iz_next)',...
                Nnext_Nxz(:),x_next_Nxzz(:,iz_next));
            S_next(:,iz_next) = interp2(N_meshgrid,x_meshgrid,S_Nxz(:,:,iz_next)',...
                Nnext_Nxz(:),x_next_Nxzz(:,iz_next));
            zeta_next(:,iz_next) = zeta_fun(Nnext_Nxz_notrunc(:), phi_next_Nxzz(:,iz_next));
            SDF_Nxzz(:,iz_next) = parms.beta*((C_next(:,iz_next)./C_Nxz(:)).^(-parms.gamma))...
                .*(zeta_next(:,iz_next)./zeta_Nxz(:));
        end
        
        CV_S(:) = sum(StateTransitionProbs_Nxzz.*SDF_Nxzz.*S_next,2);

        % check for convergence
        VFIerr = max(abs(CV_S(:) - CV_S_old(:)));
        
        idx_band = find(iter>=ITER_BANDS(1:end-1) & iter<=ITER_BANDS(2:end),1,'first');
        UPDATEWT_OLDS = UPDATEWT_BANDS(idx_band);

        CV_S_old(:,:,:) = UPDATEWT_OLDS*CV_S_old + (1 - UPDATEWT_OLDS)*CV_S;

        if VFIerr<=ERRTOL || iter>=MAXITERS
            bContinue = false;
        end
        
    end

    gesol.iterations = iter;
    gesol.VFIerr = VFIerr;
    gesol.CV = CV_S;
    gesol.S = S_Nxz;
    gesol.theta = theta_Nxz;
    gesol.C = C_Nxz;
    gesol.V = V_Nxz;
    gesol.f = f_Nxz;
    gesol.g = g_Nxz;
    gesol.spd1period = reshape(SDF_Nxzz,[NN Nx Nz Nz]);
    gesol.spd1period_format = 'N by x by z by zprime';

    % needed for bond decomposition
    SDF_Nxzz_cons = zeros(NN*Nx*Nz,Nz);
    SDF_Nxzz_zeta = zeros(NN*Nx*Nz,Nz);
    for iz_next = 1:parms.Nz
        SDF_Nxzz_cons(:,iz_next) = parms.beta*((C_next(:,iz_next)./C_Nxz(:)).^(-parms.gamma));
        SDF_Nxzz_zeta(:,iz_next) = zeta_next(:,iz_next)./zeta_Nxz(:);
    end

    gesol.spd1period_cons_comp = reshape(SDF_Nxzz_cons,[NN Nx Nz Nz]);
    gesol.spd1period_zeta_comp = reshape(SDF_Nxzz_zeta,[NN Nx Nz Nz]);
    
    gesol.wages = parms.eta*EXPZ_Nxz + (1-parms.eta)*parms.b ...
        + parms.eta*(1-parms.eta)*gesol.f.*CV_S;

end